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A central step in time-distance local helioseismology techniques is to obtain travel times of packets of wave signals 
between points or sets of points on the visible surface. Standard ways of determining group or phase travel times involve 
cross-correlating the signal between locations at the solar photosphere and determining the shift of the envelope of this 
cross correlation function, or a zero crossing, using a standard wavelet or a reference wave packet. Here a novel method is 
described which makes use of triple correlations, i.e. cross-correlating signal between three locations. By using an average 
triple correlation as reference, differential travel times can be extracted in a straightforward manner. 
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1 Introduction 

Since the first realisation of the usefulness of time-distance 
helioseismology (Duvall et al. 1993) the technique has seen 
considerable application and has developed also in terms 
of the tools used to interpret the relationship between the 
time delays measured for wave packets travelling a specified 
horizontal distance and the perturbations in temperature and 
velocity fields below the visible surface. A recent review 
describing the current state of the subject is Gizon & Birch 
2005. 

Time-distance helioseismic processing starts with the 
raw visible surface velocity data in the form of Doppler- 
grams and ends with performing the inverse problem to ob- 
tain tomographic images of sub-surface structures. An es- 
sential intermediate step is to carry out cross-correlation of 
velocity time series between locations where the relevant 
spectral line is formed, separated by one or more specified 
horizontal distances. Normally the cross-correlation is car- 
ried out between patches of pixels, e.g. centre-annulus or 
arc-to-arc cross-correlation, in order to improve the signal- 
to-noise ratio. Also, the signal is often pre-filtered, to let 
through only signals travelling at horizontal phase speeds 
confined within a small band. 

The cross-correlation function typically takes the form 
of a wavelet centered at a certain delay time away from 0, 
which is superposed on noise. Therefore one can fit e.g. a 
Morlet wavelet to this in order to extract the time of the 
maximum of the envelope, which is a proxy for the group 
travel time, or the time of a selected zero-crossing which 
is a proxy forn the phase travel time. These travel times 
are then the input for the inverse problem. Recently, an im- 
provement of this technique has been proposed (Gizon & 
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Birch 2004) where instead of a Morlet wavelet, a reference 
wavelet is used, for instance obtained from quiet Sun data, 
or computed from a reference solar model. The details of 
the fitting and time-delay extraction can be found in Gizon 
& Birch 2004. This technique is more robust in the sense 
that it continues to function at lower signal-to-noise ratios. 
Nevertheless it is still somewhat cumbersome and there is 
scope for alternative estimators for travel times. In this pa- 
per one possible alternative is discussed, which is based 
on triple correlations. In some sense the technique is car- 
ried over from for instance speckle masking interferometry 
where it has been applied succesfully for compensating for 
atmospheric seeing to achieve aperture-limited high spatial 
resolution imaging in a variety of settings (cf. Lohmann, 
Weigelt & Wirnitzer 1983, Bartelt, Lohmann & Wirnitzer 
1984). 

2 Technique 

The triple correlation of three time series /i f 2 and fa is 
defined as : 

c(n,T 2 )= f dtf 1 (t)f 2 (t + n)f 3 (t + T 2 ). (1) 

This is most conveniently calculated in the Fourier domain 
where the Fourier transform C of the correlation function 
c is related to the Fourier transforms F\ , F 2 and F3 of the 
three time series by taking their product as follows : 

C(wi,W2) = F 1 {cj 1 )F 2 {cj 2 )F;(u; 1 + lj 2 ) . (2) 

Although it is possible to calculate these triple products di- 
rectly from the Fourier transform of single pixel time se- 
ries, in the first instance the pre-processing is retained from 
the standard two-point cross correlations, using averaging 
masks to improve the signal-to-noise ratio and a phase speed 
filter to isolate a wave packet. This is done in order to com- 
pare the results of the two methods, and also in order to 
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Fig. 1 The weights in the averaging masks for the triple 
correlation shown as a linear color scale with red the small- 
est and white the largest relative weight. Note that the masks 
are all slightly displaced outward from the origin to bring 
out the structure more clearly. 



assess the sensitivity of the results to noise propagation, by 
adjusting the parameters of the masks and filters. 

A GONG (Global Oscillations Network Group) data- 
cube, tracked at the solar rotation rate, and re-sampled onto 
a common spatial grid is used for testing purposes. This 
cube is Fourier transformed in the two spatial directions and 
in the time direction. The mask averaging is done conve- 
niently in the Fourier domain, because of efficiency. If one 
needs to mask and average data using a mask with a fixed 
shape which is shifted around to sample the entire domain 
available, this can be expressed as a convolution. In the (spa- 
tial) Fourier domain the convolution becomes a simple mul- 
tiplication of the spatial Fourier transform of the data and 
the Fourier transform of the mask. In this way all necessary 
averages are obtained at once. Also, the volume of data that 
is needed for subsequent steps is reduced drastically since it 
is now significantly oversampled in the spatial directions so 
that a lot of redundant data can be dropped from the subse- 
quent steps in the analysis. 

Similarly the application of the phase speed filter to the 
data is done by Fourier transforming the phase speed filter 
and multiplying through with the Fourier transformed data- 
cube. The phase speed filter H used is Gaussian in shape : 



H{k,w) 
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With this definition A„ corresponds to the full-width at half 
maximum (FWHM). The masks chosen for the three point 
correlation function are arcs. The centres of the three arcs 
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Fig. 2 In the Fourier domain a triple correlation fills a 
hexagonal area, bounded by the Nyquist frequency in u>\ 
u>2 and u>i + u>2, shown in the figure. Since the time series 
are real there is redundancy and only two quadrants need be 
manipulated. The arrows indicate which points are complex 
conjugate pairs. 



are placed on an equilateral triangle at a distance corre- 
sponding to 29 Mm. For the purposes of these feasibility 
tests only a single distance is used. For full tomographic in- 
versions, a range of distances would be used, and the prop- 
erties of the phase speed filter and mask would be adjusted 
accordingly. The weighting is not uniform over the arcs : 
the profile of the arc in the radial direction as well as in the 
azimuth angle is Gaussian, with a FWHM of 10 Mm and 
90° respectively. The masks are illustrated in Fig. ^ which 
shows as a color scale the averaging weight for the three 
arcs. Each arc is normalised to have a unit sum over the 
pixel weights. The parameters of the phase speed filter are 
set to : 



v cen = 44 km/ s 
A„ =17.6km/s 



(4) 



Before discussing the results of the triple correlations it 
is useful to consider what one would expect. The time series 
are sampled at a finite rate, which implies that there is no in- 
formation in the Fourier domain for frequencies that are in 
absolute value above the Nyquist frequency. In the domain 
of the two Fourier frequencies u>x,u>2 (see Fig. [2) one can 
therefore restrict the analysis to a square region centered 
on the origin. Two triangles are further 'cut away' because 
these fall outside the Nyquist range for ui\ +u>2- This hexag- 
onal area contains the complex valued Fourier transform of 
the triple correlation function. Since the original time se- 
ries, and therefore also their triple correlation, is real valued, 
points in this diagram that are mirror images with respect to 
the origin (see the arrows in Fig.|2ji are complex conjugates. 



C(-Wi,-£J 2 ) = C*(wi,cj 2 ) 



(5) 



There is therefore no need to retain all 4 quadrants : two suf- 
fice and here the choice is made to retain quadrants 1 and 4. 
In speckle interferometry there are further symmetries in the 
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quantities that are being correlated, with the consequence 
that the domain necessary to fully specify the information 
can be cut down further. That is not the case here ; the time 
series are not invariant under time reversal for instance, and 
therefore no further reduction is possible. 

In the time domain, two-point correlation functions tend 
to have a wavelet shape superposed on noise. In the Fourier 
domain, two-point correlation functions are normally also 
structured with multiple peaks, and three point correlation 
functions therefore are as well. Most commonly in time- 
distance helioseismology, the shape of the wavelet is not 
used. The only parameter that is extracted is a time delay, 
either from a zero crossing or from the location of the max- 
imum of an envelope. It is consistent with this to assume 
that to lowest order the shape of the triple correlation i.e. 
its complex modulus, does not change over a field of triple 
correlations. To this order, the only quantity of interest is 
the relative displacement of this wavelet to larger or smaller 
delays, which in the Fourier domain corresponds to a shift 
in the complex phase. In other words the triple correlation 
wavelet w{t\ , t 2 ) in the time domain is presumed to be re- 
lated to an (ensemble) average (w) (ti , T2) as : 

w(n,T 2 )= J dr[ J &t' 2 (w)(t[,t 2 )x 

^-ri)^-r 2 ), (6) 

in which the 6 are Dirac delta functions. This describes a 
convolution and therefore can be expressed as a simple mul- 
tiplication in the Fourier domain. This suggests an approach 
such as Wiener filtering, in which the Fourier transform of 
the triple correlation W is divided by the Fourier transform 
of the ensemble averaged triple correlation (W). Then what 
one retains is the Fourier transforms of the two 5-functions 
i.e. a function of the form e 1 ^ in which the complex phase 
<fi is : 

4> = W1A1 + w 2 A 2 . (7) 

The Ai and A 2 are differences in travel time between and 
individual triple correlation and the average, as represented 
by the average triple correlation (w) . 

In Figs. 13 and |4] the complex modulus and phase are 
shown of the average triple correlation in the Fourier do- 
main (W) obtained by averaging the Fourier amplitudes and 
phases over all of the 8x8 triple correlations over the field 
covered by the data-cube. The phase is set to where the 
modulus is below a certain threshold. By construction the 
distances between the three points being correlated is the 
same and therefore one would expect the maximum corre- 
lation to occur for n = t 2 in the average correlation (w) . 
In the Fourier domain this corresponds to structures aligned 
along u>i + 0J2 = const.. For the same reason the individual 
cross correlations and the ratios WJ (W) should show the 
same structure. This is particularly clear in Fig. |4] In Fig. 
[3]the modulus is high only in very localised regions so that 
such a structure is lost. 

In calculating the ratio of triple correlations there is ev- 
idently a problem in large regions of the domain where the 
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Fig. 3 The modulus of the complex valued Fourier trans- 
form of the triple correlation in quadrants 1 and 4, when av- 
eraged over the entire 8x8 field of triple correlations. The 
color scale is logarithmic covering 20 powers of 10 ; black 
and red representing smallest and blue and white largest val- 
ues. 




Fig. 4 The complex phase of the Fourier transform of the 
triple correlation in quadrants 1 and 4, when averaged over 
the entire 8x8 field of triple correlations. The color scale 
is linear : yellow is zero phase, red and black is negative 
phase, blue and white is positive phase. 
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Fig. 5 The modulus of the complex valued Fourier trans- 
form of the triple correlation ratio $ in quadrants 1 and 4, 
at a location near the centre of the field. The color scale is 
logarithmic covering 5 powers of 10 with color coding as in 
Fig. 3. 

average (W) = 0. This is not an uncommon problem in in- 
versions, and one expression of the need for regularisation 
in inverse problems. In this case the most straightforward 
regularisation is similar to what is done in singular value 
decomposition. Rather than using WJ (W) one uses : 

W(UJI,LU 2 ) 



Fig. 6 The complex phase of the Fourier transform of the 
triple correlation ratio $ in quadrants 1 and 4, at a location 
near the centre of the field. The color scale is linear : yel- 
low is zero phase, red and black is negative phase, blue and 
white is positive phase. 



$(wi,w 2 ) 



(8) 



(W(ui, uj 2 )) + e 

where e is a suitably small number which needs to be ad- 
justed to the problem at hand. Here e = 10" 13 , which is 
small compared to the maximum value in the average triple 
correlation (W) max ~ 10 16 . 

If the shape of the wavelet were indeed constant over the 
field, and only displaced to larger or smaller relative time 
delays, the maps of the complex modulus of $ would be 
featureless : equal to 1 everywhere. The complex phase of 
$ would show tructures aligned along u)i + uj 2 = const, as 
mentioned above. However, the group and phase travel time 
of waves both change over the field, but not necessarily by 
the same value, due to the fact that the dispersion relation 
is complex. Only if the sub-photospheric region of the Sun 
were to satisfy : 

l k l duj 

— ^77 = const. (9) 
uj o\k\ 

would <f> show this behaviour. Since this is not the case, the 
modulus of $ also shows structure, as does its phase. 

In Figs.[5]and[6]the triple correlation ratio in the Fourier 
domain $ is shown for a single location near the centre of 
the field. Fig.|5]shows the logarithm of the modulus and Fig. 




Fig. 7 The modulus of the complex valued Fourier trans- 
form of the triple correlation ratio $ in quadrants 1 and 4, 
at the same location as Fig. 5 but without any filtering. The 
color scale is logarithmic covering 5 powers of 10 with color 
coding as in Fig. 3. 
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Fig. 8 The modulus of the complex valued Fourier trans- 
form of the triple correlation ratio $ in quadrants 1 and 4, at 
the same location as Fig. 5 but with a more localised mask. 
The color scale is logarithmic covering 5 powers of 10 with 
color coding as in Fig. 3. 



|6]shows the phase. There are two clear ridges present as ex- 
pected, both in the modulus and in the phase diagrams, from 
which a time-delay can be extracted. It is clear that along 
rdiges of uii + 0J2 — const, there is little variation other 
than what is created because of the regularisation process. 
One can therefore integrate along the ridges and reduce the 
time-delay extraction to determining a single mean value 
Ti = T2 = r m from the remaining 1-D Fourier transform. 

3 Discussion & Conclusions 

The triple correlation technique to extract relative time de- 
lays is demonstrated to function on a data-set obtained from 
GONG. The tests discussed in this paper are done on real 
data for an active region, with standard filtering and aver- 
aging applied and therefore it is evident that noise does not 
pose a significant problem in extracting time delays. How- 
ever, further work is necessary to establish the relationship 
between the relative time delays recovered in this way, and 
the perturbations in sub-photospheric layers. As has been 
pointed out in Gizon & Birch 2002 and in Jensen & Pijpers 
2003 it is possible to express this relationship in terms of a 
linear inverse problem with known kernels, as long as per- 
turbations from a known equilibrium are small. However, 
the precise form that these kernels take does depend on the 
processing in terms of averaging and filtering. Although the 
differences with existing kernels are unlikely to be substan- 
tial, appropriate kernels for the triple correlation processing 
are still to be derived. 



A point to note is that the time delays extracted are mea- 
sured relative to the mean solar sub-photosphere over the 
time period covered by the data. This mean is not necessar- 
ily identical to a standard solar model. In order to be able to 
interpret the time delays with the appropriate kernels calcu- 
lated from a standard model one should also extract the dif- 
ferential time delay of the average with respect to the stan- 
dard solar model. This can be done by calculating the triple 
correlation W^od as it would be for the model and extract- 
ing the delay from the (regularised) ratio (W) / (W moc i + e). 
The reason to proceed in this way is that the average triple 
correlation (W) is a horizontal average for the region, for 
which a comparison with a (horizontally invariant) standard 
solar model is more straightforward to interpret. The hor- 
izontal perturbations measured by individual triple corre- 
lations are more likely to be linear when compared with 
a horizontal average, than they are when compared with a 
standard solar model. 

In order to illustrate the sensitivity to the filtering of the 
triple correlations, the same data-cube was processed but 
without using any filtering at all. The equivalent of Fig. [5] 
is shown in Fig.0 It is clear that qualitatively similar struc- 
ture is present in both images so that even for unfiltered data 
it would be possible to extract a relative time delay. How- 
ever in the unfiltered processing the relative delays would 
be quite different from those recovered from filtered data. 
This is due to the fact that without the filtering there is a 
mixture of wave modes with very different depths of pene- 
tration into the sub-surface layers. The relative delays with 
respect to the mean appear to be much larger here, which 
suggests that these data are dominated by modes that remain 
in shallow layers, such as f modes. The influence of noise, 
either measurement noise or intrinsic solar noise, appears to 
be minimal. 

A separate test is to retain the filtering, but to reduce the 
extent of the averaging mask. In Fig.[8]both the radial extent 
and the azimuthal extent (FWHM) of the arcs is reduced by 
a factor of \/2. This figure shows essentially the same struc- 
ture as Fig. [5] although perhaps more of the expected peri- 
odic structure is present in the form of a third ridge, which 
implies that the time delay is somewhat better constrained 
for the smaller mask. Further optimisation of the combina- 
tion of filtering and spatial averaging is in progress. 

The computational burden of processing this data cube 
of 512 slices of 128x128 pixels is 12 min on an 8 proces- 
sor machine. This is very similar to the standard two point 
correlation processing, with a more direct and robust ac- 
cess to the time delays. The preliminary tests with filtering 
and averaging indicate that perhaps less averaging is neces- 
sary compared to standard two-point correlation functions, 
which would be beneficial for subsequent tomography. 
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